function[X] = SimulateBarristaArrivals(alpha1,alpha2,alpha3,d1,d2,T,n)
%%% Generates n arrivals in the time interval [0,T] 
%%% [0, d1]    = phase 1 (espresso), with intensity coefficient alpha1
%%% [d1, T-d2] = phase 2 (macchiato), with intensity coefficient alpha2
%%% [T-d2, T]  = phase 3 (ristretto), with intensity coefficient alpha3

%CDF until time T-d
C = (alpha1*alpha2*alpha3/T ) / ...
    ( (1-d1/T)^alpha2 *(alpha1-alpha2)*alpha3 + alpha3*alpha2*(1-d1/T)^(alpha2-alpha1) + ...
    alpha1*(alpha2-alpha3)*(d2/T)^alpha2 );

u_d1 = C*T/alpha1 * (1-d1/T)^(alpha2-alpha1)* (1- (1-d1/T)^alpha1 );
u_d2 = 1 - C*T*(d2/T)^alpha2 / alpha3;

counter1 = 0; counter2 = 0; counter3 = 0;
for i=1:n,
    u = rand(1);
    if u < u_d1,
        X(i) = T - T*( 1- u*alpha1 *(1-d1/T)^(alpha1-alpha2) /(C*T) )^(1/alpha1);
        counter1 = counter1 + 1;
    elseif u < u_d2,
        X(i) = T - T*( (1-d1/T)^alpha2  - alpha2*(u - u_d1)/(C*T) )^(1/alpha2);
        counter2 = counter2 + 1;
    else
        X(i) = T - T*( alpha3*(1-u)*(d2/T)^(alpha3-alpha2)/(C*T) )^(1/alpha3);
        counter3 = counter3 + 1;
    end;
end;

sprintf('There are %d arrivals between [0, %.2f], \n %d arrivals between [%.2f, %.5f ], \n and %d arrivals in the last interval [%.5f, %.2f]',...
    counter1, d1, counter2, d1, T-d2, counter3, T-d2, T)